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Abstract - Age-related bone loss and postmenopausal osteoporosis are disorders of bone remodelling, 
in which less bone is reformed than resorbed. Yet, this dysregulation of bone remodelling does not 
occur equally in all bone regions. Loss of bone is more pronounced near and at the endocortex, 
leading to cortical wall thinning and medullary cavity expansion, a process sometimes referred to as 
"trabecularisation" or "cancellisation". Cortical wall thinning is of primary concern in osteoporosis due 
to the strong deterioration of bone mechanical properties that it is associated with. In this paper, we 
examine the possibility that the non-uniformity of microscopic bone surface availability could explain 
the non-uniformity of bone loss in osteoporosis. We use a computational model of bone remodelling in 
which microscopic bone surface availability influences bone turnover rate and simulate the evolution 
of the bone volume fraction profile across the midshaft of a long bone. We find that bone loss is 
accelerated near the endocortical wall where the specific surface is highest. Over time, this leads to 
a substantial reduction of cortical wall thickness from the endosteum. The associated expansion of 
the medullary cavity can be made to match experimentally observed cross-sectional data from the 
Melbourne Femur Collection. Finally, we calculate the redistribution of the mechanical stresses in this 
evolving bone structure and show that mechanical load becomes critically transferred to the periosteal 
cortical bone. 

Key words: osteoporosis, endocortical bone loss, cortical thinning, specific surface, mathematical 
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1 Introduction 

It is well known that cortical porosity of bone in- 
creases in osteoporosis leading to a reduction in 
bone stiffness and strength, ultimately increasing the 
risk of fracture [1-7]. The temporal evolution of 
this deterioration process in single individuals is still 
poorly understood. Most of our current knowledge is 
deduced from cross-sectional data collected by micro- 
computed tomography of bones. These data indicate 
that changes in cortical porosity are not uniformly 
distributed across the cortical thickness. In particular, 
in long bones, loss is more pronounced near the endo- 
cortical surface [1-8]. Age-related bone loss is there- 
fore characterized by an expansion of the marrow (or 
medullary) cavity. The result of this process is often 
referred to as "trabecularisation" of cortical bone, 
since the aged cortical bone exhibits morphological 
similarities with trabecular bone [8]. In the following, 
we refer to this non-uniform bone loss of cortical 
bone as endocortical bone loss. Endocortical bone loss 
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results in a reduction in cortical wall thickness and 
in an increase in porosity which consequently reduces 
the load-bearing capacity of bone [4] . 

Several mechanisms have been hypothesised to 
drive bone loss in osteoporosis, including hormonal 
changes, reduced physical activity, and an evolving 
bone micro-structure leading to changes in bone sur- 
face availability. Hormonal changes are associated 
with increased remodelling activity, reduction in os- 
teoblast activity, reduction in osteoblast number, and 
bone imbalance within single basic multicellular units 
(BMUs) [4, 9, 10]. To explain the non-uniformity of 
age-related bone loss, RB Martin suggested that bone 
loss may be influenced by the microscopic availability 
of bone surface [11, 12]. Based on experimental 
data, Martin demonstrated a remarkable relationship 
between bone porosity and bone specific surface (i.e., 
surface area per volume of bone tissue). Using this 
relationship, Martin found from a simple mathemat- 
ical model that in cortical bone increased resorption 
activity of bone cells would lead to an increase in sur- 
face area which in turn would lead to an acceleration 
in resorption activity, hence, creating a morphological 



(or geometrical) feedback [11]. In trabecular bone 
this situation is reversed: increased resorption activity 
would lead to a decrease in surface area which in 
turn would lead to a reduction in resorption activ- 
ity. Martin's hypothesis that bone remodelling in 
osteoporosis contains such a morphological feedback 
by the evolving bone micro-structure is difficult to 
study experimentally and Martin did not provide any 
comparison to experimental data. Indeed, the time 
scales involved in osteoporosis are large, and there is 
no easy experimental control of the strength of the 
morphological feedback, nor of how bone morphology 
changes with age. 

In this paper we propose a computational mod- 
elling approach suited to the investigation of the non- 
uniformity of bone loss in osteoporosis attributed to 
a morphological regulation. Changes in cortical bone 
volume fraction occurring in the midshaft of human 
femur bone are evolved numerically assuming an 
initial bone state and an osteoporotic condition. Our 
computational model follows Martin's idea [11, 12]: 
we hypothesise that the experimentally observed non- 
uniform loss of bone in osteoporosis (e.g., more pro- 
nounced near and at the endocortex) is due to the 
superposition of: 

(i) Hormonal changes and/or changes in the overall 
loading, leading to a uniform baseline of bone 
loss across the cortical bone; 

(ii) An evolving bone microstructure, that locally 
increases or decreases bone remodelling activity 
(by morphological feedback), and so locally in- 
creases or decreases bone loss compared to the 
hormonal/mechanical baseline. 

Our main hypothesis is that the above mechanisms 
are able to explain (at least in part) the rate of en- 
docortical bone loss represented by the expansion of 
the medullary cavity, for which there are experimental 
data available [2]. The predominant loss of bone 
in the endocortical region is expected to increasingly 
transfer mechanical loading towards the periosteum. 

To test these hypotheses we extend the (purely 
temporal) mathematical model initially proposed by 
Martin [12] to include a spatial component. The 
resulting spatio-temporal description enables us to fol- 
low the evolution of the distribution of bone volume 
fraction across the midshaft of a long bone. The 
feedback of bone morphology on the bone cells is im- 
plemented using the phenomenological relationship 
between bone porosity and bone specific surface men- 
tioned above. Furthermore, we investigate changes 
in the distribution of mechanical stresses across the 



midshaft attributable to the evolving bone microstruc- 
ture. The distribution of the mechanical stresses is 
calculated using an extension of classical beam theory 
to materials of non-uniform composition. This theory 
and the mathematical model of bone remodelling are 
presented in Section 2. The numerical results are 
shown in Section 3 and discussed in Section 4. 

2 Methods 

In this section, we first present a mathematical model 
of bone remodelling well-suited to investigate the role 
of non-uniform bone surface availability for endo- 
cortical bone loss such as that which occurs leading 
to osteoporosis. We then present the method by 
which the redistribution of the internal mechanical 
stresses will be calculated when bone has evolving 
non-uniform properties. 

Computational model of bone remodelling in- 
cluding morphological regulation 

Bone remodelling is a complex metabolic process that 
involves regulations at several time scales and length 
scales. To follow the evolution of bone properties 
that are relevant to the macroscopic degradation of 
bone in osteoporosis, we consider a computational 
model focused on a millimetre-size tissue scale (see 
Figure la). This observation scale is large enough 
for local resorption and formation processes to be 
of macroscopic significance for the mechanical and 
microstructural properties of the tissue [13, 14], and 
yet small enough for the cellular origin of resorption 
and formation to be considered with biochemical and 
morphological regulations [15, 16]. At the tissue 
scale, the various cell populations involved in bone 
remodelling can be represented by continuous vari- 
ables, i.e. local cell densities n(r, t), by means of 
local averages over a so-called 'representative volume 
element' (RVE) of the tissue V T 3-8 mm 3 : 

n(r,t) = — — , (1) 

where N(r, t) is the number of cells in the volume V T 
centred at position r, at time t (see Figure la). 

Stiffness properties of a millimetre-size portion of 
bone tissue are determined to a great extent by 
the bone volume fraction (or equivalently the bone 
porosity) and the elastic properties of the mineralised 
bone matrix [13, 14]. 1 In osteoporosis, both bone 

1 Generally speaking, pore shape and orientation are other 
important factors determining stiffness properties of porous mate- 
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bone volume fraction / bm [-] 

Figure 1 - (a) Compressive force N and bending moment M acting in a cross-section (of area A) of a long bone. The representative 
volume element (RVE) of the tissue at position r serves to define local spatial averages of bone properties, such as cell densities, bone 
matrix volume fraction f bm and specific surface s v . (b) Relationship between specific surface s v and bone matrix volume fraction / bm as 
inEq. (7) (redrawn from Ref. [12]). 



volume fraction and the overall degree of mineral- 
isation of bone matrix are decreased. Osteoporosis 
is associated with increased remodelling activity with 
net bone loss per remodelling event that becomes 
progressively stronger with time [4]. This rapid 
remodelling both exacerbates bone loss and replaces 
more densely mineralised matrix with younger, less 
mineralised matrix [17]. In this paper, we focus 
on age-related changes in the morphological rather 
than mineral properties of the bone tissue. Such 
morphological changes have been studied extensively 
in cadaver bone specimens from the Melbourne Fe- 
mur Collection [2,5-7]. Changes in the degree of 
mineralisation of bone matrix will be considered in 
future studies. The bone volume fraction / bm (r, t) of a 
portion of bone tissue represents the relative amount 
of bone matrix (volume V bm ) in the RVE of volume 
V T : 2 

/bmO,t)= • (2) 

Local bone matrix is continually renewed by re- 
modelling. In homeostasis, resorption and formation 
are balanced and remodelling leads to bone turnover 
without loss or gain. 3 In osteoporosis, less bone 
is reformed than is resorbed, so bone turnover is 

rials. In cortical bone, however, pores are fairly uniformly aligned 
cylindrical Haversian canals and their spatial distribution within 
cross-sections, for example, is unimportant for tissue stiffness at 
the millimetre scale [14]. 

2 The bone volume fraction / bm is also equal to 1 — where $ 
is the bone porosity. In Ref. [15], / bm was denoted by bv. 

3 In cortical bone, remodelling generates both type I osteons, 
associated with a new Haversian canal, and type II osteons, asso- 
ciated with a pre-existing Haversian canal [18, 19]. The creation 
of a new Haversian canal in type I osteons therefore always 
implies a net bone loss [20]. However, the exact proportion of 
type I over type II osteons in normal remodelling is controversial. 
Whilst the density of vascular channels increases progressively 
with age [19], age-related bone loss is associated with increased 



associated with a net bone loss. The evolution of the 
local bone volume fraction _f bm is governed by (i) the 
local densities of active osteoblasts (OB a ) and active 
osteoclasts (OC a ) and (ii) these cells' activity rates: 

j- t fbm(r, = fcformOBaO, t) - fc res OC a (r, t), (3) 

where fc form is the volume of new bone synthesised 
per unit time by a single active osteoblast and fc res 
the volume of bone resorbed per unit time by a single 
active osteoclast. 

Equation (3) expresses generally the balance be- 
tween bone formation and bone resorption: fcf or m OB a 
corresponds to the formation rate (bone volume frac- 
tion synthesised per unit time) and ?c res OC a corre- 
sponds to the resorption rate (bone volume fraction 
resorbed per unit time). The difficulty in Eq. (3) lies in 
determining the evolution of the populations of active 
osteoclasts and active osteoblasts, and in particular, 
in determining how their densities depend themselves 
on the presence of bone matrix. Osteoclasts and 
osteoblasts usually require a pre-existing bone sub- 
strate to conduct their resorbing and synthesising 
activities. 4 The densities of osteoclasts and osteoblasts 
thus depend on the local amount of bone surface (area 
S) available to them in the RVE, i.e., on the specific 



pore area rather than increased pore density [7]. The importance 
of cortical bone loss with age due to normal remodelling is thus 
difficult to estimate. In the present model, we do not take this 
effect into account and we will assume the same baseline of bone 
imbalance per remodelling event in osteoporosis in cortical and 
trabecular bone for simplicity. 

4 An exception is the growth of so-called 'membrane bone', 
which forms de novo without cartilaginous substrate or scaffold 
in the flat plates of the skull. 
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surface s v (r, t): 5 

s v (r,t) = ——. (4) 

\/rj. 

Comprehensive cell population models of os- 
teoblasts and osteoclasts including biochemical cou- 
pling have been proposed in the literature [15, 21- 
25]. In Ref. [25], the potential influence of the 
specific surface s v on various stages of osteoblast and 
osteoclast developments is investigated in the bone 
cell population model of Refs [15,24]. Here, how- 
ever, we will consider a simpler model of bone cells, 
originally formulated by Martin [12]. The advantage 
of this simpler model is that the dependence upon 
s v of the densities of active osteoclasts and active 
osteoblasts appears explicitly. Whilst full biochemical 
coupling between cells is not explicit in this model, 
the non-uniformity and evolution of the bone surface 
availability are accounted for, which are important 
to capture non-uniform effects of remodelling during 
age-related bone loss according to our hypotheses. 

Bone tissues with a large specific surface exhibit 
more remodelling, and so contain more active os- 
teoblasts and active osteoblasts, than tissues with a 
small specific surface (Figure la) [12]. Denoting by 
A OBa and A 0Ca the fractions of the bone surface at 
which there is osteoblastic and osteoclastic activity 
due to remodelling, and by cr 0Ba and <r OCa the surface 
densities of active osteoblasts and active osteoclasts 
at remodelling sites, the densities of active osteoblasts 
and active osteoclasts are given by OB a = A 0B <r 0B s v 
and OC a = A OCa cr OCa Sy. Substituting these expressions 
in Eq. (3), one thus has [12]: 

^/bm( r » — (^form^OB^OBa _ ^res^oc^ocj 5 v( r » 0- 

(5) 

Whilst Eq. (5) holds generally, in the following, it 
is assumed that fc form , ?c res , A ov A 0Ca , cr OBa , cr 0Ca 
can be taken uniform and constant in time. This is 
clearly a simplification, as the microscopic availability 
of bone surface may influence the recruitment and 
development of bone cells in a nontrivial way. These 
recruitment and development processes determine the 
fractions A OB and A oc , which would thus normally be 
expected to depend additionally on s v [25]. Here, a 
simple extensivity of cell densities in s v is assumed 
when taking OB a and OC a just proportional to s v . 

In osteoporosis, a net bone loss is associated with 
each remodelling event. To retrieve physiological 
overall rates of bone loss, Martin assumed a small 

5 We do not include surfaces of the canalicular system in this 
definition of the specific surface s v . Canalicular surfaces are not 
available to osteoclasts and osteoblasts for remodelling. 



imbalance between formation and resorption, set- 
ting [12]: 

fcform^oB a ^oB a - fcres^oc a °-oc a = -2 jim/year. (6) 

This constant baseline of bone loss means that a 2 \im- 
thick layer of bone matrix is resorbed each year on all 
bone surfaces. This baseline of bone loss is weighted 
by the specific surface to correspond to volumetric 
losses in the local RVE of the tissue V x . 

Bone tissues exhibit a wide range of microstructures 
each characterising a particular volume fraction and 
specific surface. Cortical bone typically has volume 
fractions / bm fa 0.85-0.95 with pores mainly consti- 
tuted of Haversian canals. Trabecular bone typically 
has volume fractions _f bm sa 0.15-0.55 with 'pores' cor- 
responding to the marrow space around the trabecular 
plates and struts. Based on measurements performed 
in a large variety of bone tissues, Martin [12] has 
proposed that bone volume fraction _f bm and specific 
surface s v follow an 'intrinsic' or 'universal' relation- 
ship, approximated by the following polynomial: 

s v(/bm) — a l/bm + a 2/bm 2 + a 3/bm 3 + a 4-/W 4 + a 5.fbm 5 > 

(7) 

with 

a x = 14.1/mm, a 2 = — 10.5/mm, a 3 = — 17.8/mm, 
a 4 = 43.0/mm, a 5 = -28.8/mm. (8) 

The relation (7) is plotted in Figure lb together with 
experimental data assembled in Ref. [12] from various 
types of human bones (femur, iliac crest, vertebra, rib) 
both in health and disease. 6 Importantly, the specific 
surface exhibits a maximum s v * 4.2/mm at a bone 
volume fraction _f bm * 0.63, intermediate between 
cortical and trabecular bone. 

Using the phenomenological relation (7) and 
Eq. (6) in Eq. (5), one is able to calculate the evolution 
of bone volume fraction / bm (r, i n every region 
r of the bone tissue, given an initial distribution 
of bone volume fraction / bm (r, to)- This evolution 
enables us to determine how the medullary cavity 
expands in time for comparison with experimental 
data [2] . It is well-known that the distinction between 
cortical and trabecular bone is not straightforward in 
osteoporosis [8]. Similarly, measuring the extent of 
the medullary cavity relies on a visual or morpho- 
logical threshold to define the endocortical surface. 
For the purpose of this paper, we assume that the 

6 This figure is redrawn from Ref. [12] with the bone volume 
fraction / bm in lieu of the porosity 1— f hm . The coefficients a lt a 5 
of the phenomenological polynomial have been slightly adjusted 
from [12] such that s v (0) = s v (l) = 0. 
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medullary cavity is defined by the region of bone 
tissue in which bone volume fraction is lower than 
a threshold f^ m *. We will take this threshold as the 
bone volume fraction at which the specific surface is 
maximum according to Eq. (7), i.e. / bm * 0.63. The 
effect of choosing different threshold values for the 
determination of medullary cavity radius will also be 
investigated (see Figure 2d). 

The simplicity of the model considered here (as 
opposed to more detailed models of bone cell devel- 
opments such as that of Ref. [25]) has two practical 
consequences: (i) the specific surface s v enters as an 
explicit dependence in the right hand side of Eq. (5). 
This allows the formulation of a semi-analytical so- 
lution for the evolution of bone volume fraction f bm 
(see Appendix A); (ii) the governing equation for the 
bone volume fraction is self-contained: the state of 
the system depends only on the current distribution 
profile of the bone volume fraction. This enables us 
to easily regress the system backwards in time and to 
deduce the past distribution of bone volume fraction 
that has led to an experimentally observed medullary 
cavity expansion (see Appendix B). 

Internal mechanical stress distribution 

Mechanical loading is carried non-uniformly by bone 
tissues. In particular, the cortex carries the majority 
of the loads applied to a bone [8]. Here, we consider 
a long bone subject to a compressive resultant force 
JV(r) = — IV(t)x along x and a resultant bending mo- 
ment M(t) = M y (t)y + M z (t)z in the cross-sectional 
y = (y,z) plane (see Figure la; x, y, and z denote 
units vectors along the x, y, and z axes). To estimate 
how internal mechanical stresses are redistributing 
across the midshaft section of a long bone during 
age-related endocortical bone loss, we use a simple 
theory of elasticity extended to beams of non-uniform 
composition (see, e.g., Refs [26-28]). This theory is 
based on two main assumptions: 

1 . Hooke's law is valid locally at the tissue scale, i.e., 
the stress tensor cr(r, t) and strain tensor e{r, t) 
are related by a local stiffness tensor £(r, t) of 
the tissue: 



a{r,t) = E{r,t)e{,r,t) 



(9) 



2. Near the midshaft of a long bone, small defor- 
mations generated by compression and bending 
keep the initial cross-sections planar and normal 
to the neutral axis. This geometrical constraint 
on the deformations is called the Euler-Bernoulli 
kinematic hypothesis. 



From a structural point of view, long bones near the 
midshaft can be regarded as weakly deformed long 
beams. In absence of torsional loads or twisting 
along the beam axis, as assumed here, the defor- 
mation state of such beams is well approximated by 
the Euler-Bernoulli hypothesis [26-28]. The Euler- 
Bernoulli hypothesis implies that the strain tensor 
e(r, t) reduces to the single nonzero scalar compo- 
nent £ xx (r, t) [27,28]. For an orthotropic material 
such as bone [13], this implies by Eq. (9) that there 
are no shear stresses. The only nonzero components 
of the stress tensor for such materials are the normal 
stresses cr xx , cr yy and ct zz . The normal stresses 
(jyy and ct zz are induced by the 'Poisson effect' (i.e., 
thickening of a material under compression) and are 
usually small except in particular materials such as 
rubber [27,28], [29, §17]. The normal stresses a yy 
and a zz do not participate directly to the transfer 
of the resultant compressive force iV(t) and bending 
moment M(t) across the tissue. Indeed, the resultant 
force iV(t) and moment M(t) are given as integrals 
of the stress distribution a xx {r, t) in the midshaft 
cross-section [27, Sees 5.3, 6.2] (see below). For this 
reason, in the following we will focus on the distribu- 
tion of normal stresses <j xx {r, t) only. From Hooke's 
law (9), one has a xx {r, t) = E mi (r , t)e xx (r , t), and 
so only the single component E nu of the stiffness ten- 
sor needs to be considered, which corresponds to the 
compressive stiffness. 7 We will omit any indices for 
notational simplicity, and write £ ull (r, t) = E(r, t). 

We consider a cross-section at a position x along the 
bone axis, near the midshaft, and use the coordinates 
y = (y,z) to denote a position in the cross-sectional 
plane (see Figure la). By definition, the resultant 
force iV(t) and resultant moment M(t) are given as 
integrals of the stress distribution cr(y, t) = o xx (r, t) 
in the midshaft cross-section (we omit both the cross- 
section position 'x' and the component indices 'xx' 
from the notation). With Hooke's law (9), one thus 
has [27, Sees 5.3, 6.2]: 



E{y,t)e{y,t)c\yc\z, 



(10) 



J A 



MJt)=- 



zE{y,t)e(y,t)dydz, (11) 



M z (t) : 



yE{y, t)e{y,t)dydz, 



(12) 



where the integrals are carried over the bone cross- 
sectional area A (Figure la). The geometric constraint 



7 The calculation of the other normal stresses a yy and a zz 
depends on other components of the stiffness tensor. Using bone's 
orthotropic stiffness property [13] and the fact that only e xx =fc 0: 
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imposed by the Euler-Bernoulli condition implies that 
the strain distribution e{y, t) is a linear function of 
the transverse spatial coordinates y = (y,z) (see e.g. 
Ref. [27, Sec. 5.2]): 

s(y,t) = s 1 (t)-K 3 (t)y + K 2 (t)z, (13) 

where e 1 is the sectional axial strain, and k 2 and k 3 
are the sectional beam curvature about the y and z 
axes, respectively. At each time t, the Euler-Bernoulli 
condition (13) together with the constraints (10)- 
(12) form a system of three linear equations deter- 
mining the three unknowns £i(t), K" 2 (t) and fc 3 (t). 
The coefficients of this system of equations involve so- 
called 'sectional stiffness coefficients', i.e., zero order, 
first order and second order moments of the non- 
uniform stiffness across the section [27, Sec. 6.2]: 



the onset of osteoporosis. The time elapsed since 
onset of osteoporosis is denoted by At p = t — 1 . 
A constant compressive force iV(t) = 0.4 kN and 
constant bending moment around the z axis M z (t) = 
2 kN • m, M y (t) = 0, are assumed when calculating 
the redistribution of the mechanical stresses. Below, 
we first show how the bone volume fraction /bm(P> 
evolves by means of Eq. (5) from a typically expected 
radial profile, and we determine the medullary cavity 
expansion that this evolution induces (Section 3.1). 
In a second numerical experiment, we calculate the 
"inverse problem", i.e., we use morphological data 
by Feik etal. [2] compiled into an experimentally 
observed average medullary cavity expansion in men 
to determine what initial bone state would have gen- 
erated this medullary expansion (Section 3.2). 



S(t): 



E(y,t) dydz, Sj(t) 



Ja 



x ; £(y,t) dydz, 



I jk = \ Xj x k E(y,t)dydz, i,j = 2,3 (14) 

I A 

(with x 2 = y and x 3 = z). Once ^(t), fc 2 (t) and 
K 3 (t) are known, the mechanical stress distribution in 
the cross-section is given by use of Hooke's law: 

a(y, t) = E(y, t) [ £l (t) - K 3 (t)y + /c 2 (t)z] (15) 

The sectional stiffness coefficients (14) entering the 
system of equations to solve are specified at each time 
t provided the stiffness coefficient E(y, t) is known. 
Measurements of the compressive stiffness of bone in 
relation to the amount of bone matrix suggest that 
stiffness and bone volume fraction are related by a 
nonlinear phenomenological relationship [20, 30-32, 
and Refs cited therein], which we take here as: 

E(y, t) = £(/ bm (y, t)) = C f hm (y, t) 3 , (16) 

where C *w 15 GPa and _fbm( v > t) is governed by the 
balance equation (5). 

3 Results 

Equation (5) determines the evolution of the bone 
volume fraction /bm( r > at each position r within 
bone from a given initial condition / bm (r, t ). Here for 
simplicity only rotation-symmetric initial distributions 
are considered in a cross-section at x, implying that 
/b m stays rotation-symmetric at all times: /bm( r > = 

AmCy.z.O = /bm(p.t). where p = \/y 2 +z 2 is 
the radial coordinate. The initial condition at time 
t is assumed to correspond to the bone state at 



3.1 Evolution of a typical bone volume frac- 
tion distribution 

In Figure 2a, the bone volume fraction radial profile 
/bm(P> t) i s shown every 10 years from the onset of 
osteoporosis (At 0P = 0) for 50 years. The assumed 
initial bone volume fraction profile at At op = (solid 
line) consists of (i) a sigmoid function of p with value 
zero for p < 3.3 mm (medullary cavity), sharply 
increasing (endocortical region) to the value 0.75 at 
p 6.3 mm and (ii) a linear function more slowly in- 
creasing (intracortical region) to the value 0.99 from 
p Pa 6.3 mm to p = 12 mm (periosteum, bone radius). 
Figure 2b and 2c show the corresponding changes 
in the radial profile of internal stresses o(p,6,t) 
and specific surface %(p, t) across the midshaft. In 
contrast to / bm and s v , the internal stresses a are not 
rotation-symmetric due to the action of the bending 
loading condition. The evolution of u is shown at 
three angles 9 of the polar coordinate system (6 = 
0,^, and 7i), where 6 is measured from the y axis 
(see Figure la). The full polar distribution of stresses 
is shown in Figure 2e at At 0P = 50 years. 

The bone volume fraction in Figure 2a is seen to 
decrease in all regions of the bone cross-section due 
to the simulated osteoporotic condition. However, as 
expected from Eq. (5), bone loss is more pronounced 
in the endocortical-to-intracortical region, where the 
specific surface is high. High bone volume fractions 
near the periosteum are comparatively more pre- 
served. In the endocortical region, the bone volume 
fraction profile is both lowered and shifted towards 
the periosteum by the progression of osteoporosis, 
resulting overall in a more gradual increase of bone 
volume fraction from the medullary region to the 
periosteum than the initial profile. 
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Figure 2 - (a)-(c) Calculated radial profiles shown every 10 years since onset of osteoporosis for (a) bone volume fraction; (b) internal 
stresses at three different polar angles; (c) specific surface. The arrows indicate the evolution of the profiles with osteoporosis; (d) 
Calculated expansion of the medullary cavity radius. The effect of the threshold used to defined the radius is shown. Thresholds from 
/ bm = 0.1 to / bm = 0.9 are used (thin solid lines) in addition to the threshold f bm = / b * m at which s v is maximum (thick solid line). 
The filled region is indicative of the reduction in cortical wall thickness; (e) Stress distribution across the midshaft cross-section for 
compression and bending, exhibiting a neutral axis at y « —10 mm; (f) Contour plot of the difference in internal stresses between 
intial bone state and bone state reached after 50 years of osteoporosis. Red tones: decrease in stresses; Blue tones: increase in stresses; 
Thick grey curve: no change in stresses (zero isoline). 



7 



The radial profile of the specific surface (Figure 2c) 
gradually evolves such that its maximum originally 
at p ^ 5.5 mm is shifted towards the periosteum. 
The position of this maximum corresponds to the 
intersection between the bone volume fraction and 
the constant value f bm * 0.63 in Figure 2a because 
f bm * maximises the function s v (f bm ) in Eq. (7). Taking 
f bm * as the threshold value to define the transition 
between medullary cavity and cortex, this shift of the 
maximum of the specific surface towards the perios- 
teum corresponds to an expansion of the medullary 
cavity radius with age. This expansion is plotted 
separately in Figure 2d (thick solid line). Thin solid 
lines in Figure 2d correspond to the expansion of a 
medullary radius defined by other threshold values of 
_f bm as indicated by the labels. In Figure 2d, the onset 
of osteoporosis was assumed at 40 years of age. For 
the threshold value _f bm *, the increase in medullary 
radius is initially slow, then accelerates at age ^ 55. 
The cortical wall thickness corresponds to the distance 
between periosteal bone radius (here assumed fixed at 
12 mm) and medullary cavity radius, and is shown by 
the grey-shaded area. The expansion of the medullary 
cavity is therefore associated to a thinnning of the 
cortical wall thickness. 

These changes in morphological parameters are 
associated with a gradual transfer of the internal 
mechanical stresses towards the periosteum for = 
and 6 — n/2 (Figure 2b), whilst the stress radial 
profile at = 7i changes little. Interestingly, Figure 2b 
suggests that for each polar angle 6, there is a well- 
defined radius below which stresses are decreased 
and above which stresses are increased. The polar 
dependence of this radius can be appreciated in Fig- 
ure 2f. In Figure 2f, isolines of the difference between 
initial stress distribution and stress distribution after 
50 years of osteoporosis are shown. The radius 
below which stresses are decreased and above which 
stresses are increased corresponds to the intersection 
between <r(p, 0, At 0P = 0) and cr(p, 9, At 0P = 50) in 
Figure 2b, and to the zero isoline in Figure 2f, shown 
as a thick grey line. The region where stresses are 
increased is shaded with blue tones and the region 
where stresses are decreased is shaded with yellow- 
red tones. 

3.2 Determination of past bone volume frac- 
tion profile from an observed medullary 
cavity expansion 

As explained in Appendix B, it is possible to use our 
model to regress the system backwards in time. This 
can be used to take as input a given medullary cavity 



expansion and deduce from it the past distribution of 
bone volume fraction that leads to such an expansion 
when evolved by the model. Mathematically, this 
"inverse problem" can only be solved provided the 
medullary cavity radius is monotonously increasing. 
Also, the bone volume fraction profile can only be 
determined between the minimum and the maximum 
medullary cavity radii of the given expansion. In 
Figure 3b, we show an observation of medullary 
cavity expansion obtained by Feik et al. [2] from cross- 
sectional samples of the Melbourne Femur Collection. 
In the original data, samples were grouped in decade 
age categories (such as 31-40 years, 41-50 years etc.) 
and the medullary area was measured on each sam- 
ple. Here we have transformed medullary area into 
medullary radius by assuming a rotation-symmetric 
medullary cavity and we have interpreted the value 
in a decade age category as the value at a single 
age in the timeline (namely, the middle age of the 
category: 31-40 years — > 35 years). Because the 
"inverse problem" requires a monotonously increasing 
radius, we took the data presented in Figure 5 of 
Ref. [2] for males and discarded the first data point, 
as indicated by the hatched area in Figure 3b. 

In Figure 3a, we show the bone volume fraction 
radial profile at 35 years of age (solid line) that 
reproduces the medullary cavity radius of Figure 3b 
when evolved to 95 years (double dashed line). As in 
Figure 2a, the calculated bone volume fraction profile 
is shown every 10 years and the medullary cavity 
radius is determined by the intersection of the bone 
volume fraction profile with the constant value / bm * 
(dotted horizontal line). In Figure 3b, the minimum 
radius is ss 6.4 mm and the maximum radius is & 
9.6 mm, and so the bone volume fraction profile in 
Figure 3a is only determined between these radii as 
indicated by the excluded hatched areas. 

4 Discussion 

The study of age-related endocortical bone loss in 
humans is a challenging problem, with a multitude of 
possible causes, including biochemical, biomechanical 
and morphological regulations of bone remodelling. 
Whilst micro-CT imaging technology can give detailed 
insights into the bone microstructure at a given mo- 
ment, to follow the evolution of bone microstructure 
across a whole cross-section with a degree of detail 
and a frequency that are sufficiently informative (i) 
is not currently possible without harmful radiation 
and (ii) would require several decades of study. Fur- 
thermore, soft tissues and bone cells are not seen 
in micro-CT bone scans. Two different approaches 
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Figure 3 - (a) Determination of bone volume fraction profiles in the endocortical region (non-hatched area) from a given expansion 
of medullary cavity radius (b). The data in (b) is taken from male bone samples in Ref. [2, Fig. 5]. (The parallel periosteal expansion 
found in these male bone samples is not considered; a constant periosteal bone radius of 12 mm is assumed here instead.) Only the 
monotonously increasing part of the curve (non-hatched area) can be used to determine the evolution of bone volume fraction profiles 
(see Appendix B) . 



have been undertaken by biologists. Animals models 
of osteoporosis have been proposed, for which both 
micro-CT scans and histomorphometric analyses of 
bone samples can be performed, providing both bone 
morphological parameters and biochemical informa- 
tion. However, particularly in small animals that lack 
secondary osteons (such as mice and rats), the patho- 
genesis of osteoporosis may be quite different from 
that in humans. In humans, loss of endocortical bone 
in osteoporosis is believed to progress by expansion 
and coalescence of Haversian pores near the endos- 
teum [5-7]. A second approach is to collect cross- 
sectional data from bone samples (either ex vivo or 
post mortem) and reconstruct pseudo time sequences 
by grouping the data by age categories [2,3,5-7]. The 
advantage of this approach is the direct measurement 
of human morphological properties of bone. The 
disadvantage is that measurements can often only be 
done at a single point in time in an individual. 

To our knowledge, there is currently no experi- 
mental work specifically studying the influence of the 
microscopic bone morphology on the development 
and/or activity of bone cells. 8 Such a morphologi- 
cal feedback on bone remodelling is challenging to 
investigate experimentally for the following reasons: 
(i) a morphological feedback cannot be inhibited 
to estimate its influence (in contrast to biochemical 
regulations that can be selectively inhibited by gene 
knock-outs in animal models); (ii) there is no easy 



8 Some effects of the curvature of the bone surface on os- 
teoblast activation and/or osteoblast activity are known to occur. 
In vivo, the refilling rate of cortical bmus depends on the current 
radius of the closing cavity as evidence by double tetracycling 
labelling techniques [33-35]. In an in-vitro setup, the local 
curvature of the bone substrate was also shown to influence the 
rate of bone formation [36] . 



control of the bone microstructure in vivo; (iii) in 
osteoporosis, changes in bone volume fraction occur 
over long time scales. Computational modelling is a 
powerful tool that enables the testing of hypotheses 
for the underlying mechanisms responsible for age- 
related endocortical bone loss. Most importantly, com- 
putational modelling enables the extrapolation of data 
measured at a single point in time into a predicted 
evolution over several decades. This is of particular 
importance for the development of patient-specific 
assessment and treatment tools for osteoporosis. 

The main result of this paper is the observation that 
a morphological feedback of the bone microstructure 
(specifically, of the bone specific surface) on the bone 
cells is able to explain (at least partly) the predomi- 
nant loss of bone in the endocortical region. The non- 
uniformity of this loss can develop from a uniform 
baseline of bone imbalance per remodelling event 
as caused by systemic hormonal changes or overall 
changes in mechanical loading due to reduced physi- 
cal activity. The loss of endocortical bone leads to an 
expansion of the medullary cavity that is compatible 
with experimentally observed rates. Furthermore, the 
deterioration of bone in the endocortical region leads 
to a transfer of the mechanical stresses towards the 
periosteum. 

The evolution of the radial profile of bone volume 
fraction in the bone midshaft in Figure 2a is gov- 
erned by a partial differential equation in time only. 
However, spatial non-uniformities enter the equation 
by means of the local availability of bone surface 
embodied by the specific surface s v . Age-related 
cortical bone loss results in both an overall decrease 
in cortical bone volume fraction, and a decrease in 
cortical wall thickness due to pronounced bone loss at 



9 



the endocortical surface, where the specific surface is 
high. This is consistent with the experimental findings 
ofRef. [3]. 

The precise form of the expansion of the medullary 
cavity radius shown in Figure 2d depends strongly on 
the assumed initial bone volume fraction profile and 
on the threshold value of bone volume fraction used 
to define the boundary between medullary cavity and 
cortex. The initial slow increase in medullary radius 
in Figure 2d is due to the loss of bone in the region 
3.3 mm < p < 6 mm (initial endocortical region). In 
this region, the intial bone volume fraction increases 
sharply from to 0.75 with the radial coordinate 
p, forming a "step" (Figure 2a). Once the sloped 
"plateau" of bone volume fraction for p > 6 mm 
(initial intracortical region) has reached the threshold 
bone volume fraction defining the medullary cavity 
boundary, the expansion of medullary radius increases 
rapidly in Figure 2d before gradually slowing down 
due to the relative preservation of bone volume frac- 
tion near the periosteum. In fact, choosing a dif- 
ferent threshold volume fraction essentially delays or 
advances this general behaviour. From these results, 
one can deduce that for the influence of a morpho- 
logical feedback, the sharper the transition between 
medullary cavity (zero bone volume fraction) and 
cortex (high bone volume fraction), the slower the 
medullary cavity expansion. This observation em- 
phasises the importance of building up "good bones" 
during growth, with high intracortical volume frac- 
tions [37]. 

We note that an acceleration of the medullary cavity 
expansion is seen in males between the age groups 
81-90 and 91-100 in the data reported in Ref. [2, 
Fig. 5b] (corresponding to the age interval 85-95 
in Figure 3b). This acceleration could be explained 
by the specific shape of the bone volume fraction 
radial profile at the onset of osteoporosis. From our 
model, the bone volume fraction at age 35 shown 
in Figure 3a, calculated such that it leads to the 
experimental medullary expansion of Figure 3b, ex- 
hibits a similar sharp "step" (associated with a slow 
medullary expansion) followed by a "sloped plateau" 
(associated with a fast medullary expansion) as in 
Figure 2a. However, other physiological factors not 
accounted for in our model may play a role in this 
experimentally observed acceleration of the medullary 
expansion, such as a difference in hormonal imbal- 
ance with advancing age that could lead to an in- 
creased bone imbalance per remodelling event. We 
also note that in Figure 5 of Ref. [2], this acceleration 
of the medullary expansion occurs in parallel to an 
expansion of the subperiosteal cross-sectional area. 



Cortical wall thickness decreased in the male bone 
samples of Ref. [2] between the age groups 61-70 
and 81-90 (corresponding to the age interval 65-85 
in Figure 3b) due to both an increase in medullary 
cavity and decrease in subperiosteal cross-sectional 
area. However, cortical wall thickness was relatively 
preserved between the age group 81-90 and 91-100 
(corresponding to the age interval 85-95 in Figure 3b) 
as a result of periosteal apposition. In our model, 
changes in periosteal bone radius are not considered. 
Whilst bone remodelling processes are observed at the 
periosteum [4,38], periosteal apposition is thought 
to occur via bone modelling processes, which are not 
well described by our model in its current form. 

Figures 2b and 2e clearly show that internal stresses 
become increasingly and disproportionately redis- 
tributed on the strongest part of the bone midshaft, 
i.e., towards the bone periphery. This dispropor- 
tionate redistribution could put the bone integrity at 
increased risk. Indeed, the periosteum is a surface 
with low remodelling rate [4,38]. An acceleration 
of microcrack generation at the periosteum can be 
expected, increasing the risk of a macroscopic frac- 
ture. However, this increased risk may be offset by 
periosteal bone apposition [2,4,39]. 

Increased stresses towards the periosteum may also 
account for a conservation of intracortical pore area. 
In the current model a mechanical feedback to form 
bone at the periosteum and to slow down bone 
imbalance at the endocortex is not explicitly taken 
into account. 9 Figure 2f suggests that a potential 
mechanical feedback based on internal stresses would 
lead to further non-uniformities in the evolution of the 
distribution of bone volume fraction in the midshaft. 
Indeed, the region in which stresses are increased 
and the region in which stresses are decreased by 
the simulated osteoporotic condition have a rotation- 
asymmetric boundary as shown by the thick grey 
line in Figure 2f. The mechanical feedback would 
be particulary weak near the mechanical neutral axis 
induced by the bending loading condition and seen in 
Figure 2e near y ^ — 10 mm. Figure 2f also suggests 
that bone loss would be accelerated by a mechanical 
feedback in the endocortical region where stresses 
are decreased, concurrently to an induction of bone 
formation at the periosteum in this area, suggesting 
an outward shift of the cortex and so an evolving bone 
shape. 

In this work, bone tissue was regarded as being 



'Setting the bone imbalance to —2 u,m/year in Eq. (6) could 
be assumed to implicitly account for a mechanical feedback at 
the endocortex, but such an imbalance cannot account for bone 
formation at the periosteum. 
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solely composed of pores (vascular porosity) and 
mineralised bone matrix. In particular, the volume 
fraction of osteoid was not considered. Osteoporosis 
is an evolving bone disorder, that transitions from 
high turnover rate with low imbalance to mid-to- 
high turnover rate (possibly induced by secondary 
hyperparathyroidism) with a gradually increasing lack 
of formation [4,9]. These differences in turnover 
rates modify the local volume fraction of osteoid. 
Increased turnover rate temporarily lowers the degree 
of mineralisation of bone, which also contributes to a 
reduction in the stiffness properties of bone [13, 17]. 
Here, only the morphological dimension of the pro- 
gression of osteoporosis was considered. We also 
note that osteoid makes asymmetric the availability 
of bone surface between osteoblasts and osteoclasts. 
Indeed, active osteoclasts rarely resorb bone covered 
by osteoid, but active osteoblasts deposit osteoid on a 
previously laid osteoid layer. 

In conclusion, we found that the microscopic avail- 
ability of bone surface within the tissue, needed for 
resorption and formation processes to occur, is po- 
tentially a significant factor in focal bone loss pre- 
dominantly occuring near the endocortical wall in 
osteoporosis. Furthermore, the gradual redistribu- 
tion of the internal stresses towards the periosteum 
is suggestive of the possibility to induce modelling 
responses at the periosteum by mechanical feedback. 
Age-related endocortical bone loss involves a variety 
of biochemical and biomechanical processes that were 
not explicitly included in the model. Our model en- 
ables the extraction of the influence of the microscopic 
availability of bone surface. This influence is not 
easily assessed experimentally. Potential differences 
between experimental data on the evolution of bone 
volume fraction and the results of the model can 
therefore be attributed to these other influences, in 
particular a biomechanical feedback inducing other 
non-uniformitites in the evolution of bone volume 
fraction. 
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Appendix A Semi-analytic implicit 
solution 

The governing equation of the bone volume fraction 
has the form: 



g- t fbm( r > = -V SvC/bmO, 0)- 



(17) 



with rj = 2pim/year. This equation can be integrated 
in time at each location r : 



r 1 



dt 



(g/gQ/bm 
s v(fbm) 



f /bmO.t) 



/bmO,t ) 



df 

svtf) 



= -Tj(t-t ). 



(18) 



Denoting by g(f) an anti-derivative of l/sy(/), one 
can thus write the solution in the following implicit 
form: 

g C/bmO, 0) - g (/ bm (r , t ) = -tj (t - t ). (19) 

When s y (/) is a polynomial, one may use the repre- 
sentation: 



i " 

— =y 



A; 



(20) 



where A ; are (real) constants and f t the (possibly 
complex) roots of s v (f), with f 1 = and f 2 = 1. 
Depending on the degree of the polynomial fit used 
for Sy(/), these roots can be determined either nu- 
merically or analytically. For the particular fifth order 
polynomial proposed in Eq. (7) one has (with A; in 
millimetre) : 

fi = 0, f 2 = l, -0.611, 

/ 4 sa -0.55 -0.7 i, f 5 « 0.55 + 0.7 i, 
A x sa -2.04, A 2 «s 0.89, A 3 w 0.55, 

A 4 ** 0.30- 0.63 i, A 5 « 0.30 + 0.63 i (21) 

For polynomial s v (/), the function g(f) has thus the 
general form: 

n n 

g(f) = A; hiCf -ft) = In ( - f^), (22) 



i=l 



i=l 



and the bone volume fraction /bm( r > can be de- 
termined by solving the implicit algebraic equation 
(using Eq. (22) in Eq. (19)): 



II (Mr, t)~fi) Ai = II C/bmCr, t ) -/ ; ) Ai e"" (t -^ 

(23) 



i=l 



i=l 
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This equation can also provide an implicit equation 
for the evolution of the medullary cavity. Assuming 
a cylindrical geometry, the medullary cavity radius 
p*(t) at time t is defined by /b m (P*(' : )> = f*> 
where / * is the value of bone volume fraction chosen 
to define the endocortical wall, for example, /* = 
Am* 0-63 at which the specific surface is maximum. 
Evaluating Eq. (23) at p = p*(t), one thus has: 

n ifup%tio)-fi) Ai = n cr -/o Ai e*^ 

i=l i=l 

(24) 

Appendix B Determination of the en- 
docortical volume frac- 
tion profile at an earlier 
age 

Reversing time in the governing equation for bone 
volume fraction (17) provides the evolution equation 

f t fUr,i) = +r]s v (fUr,i)), (25) 

i.e., bone volume fraction increases with reversed time 
i — — t. Because f bm is the only state variable in 
the model presented here, one may determine the 
history of bone volume fraction from its current profile 
only. When other variables participate in determining 
the state of the system (such as bone cell densitites 
in more comprehensive models of bone cell popula- 
tions), this reversal is more complex and also likely to 
be numerically unstable. 

Below, we show that it is in fact sufficient to know 
the expansion of the medullary cavity in time to 
deduce the bone volume fraction profile at an earlier 
age (in particular, at the onset of osteoporosis) in 
the endocortical region. (This property is used in 
Figure 3 with the evolution of the medullary cavity 
area estimated in Ref. [2, Fig. 5].) For simplicity, a 
rotation-symmetric cylindrical bone shaft is assumed. 
From Eq. (24), it is clear that if the medullary cavity 
radius p*(t) is a monotonic increasing function of 
t > with values in the range [po,Pi], one can solve 
Eq. (24) for / bm (p,0) for all p e [p ,Pi]- Indeed, 
it is sufficient to find the (unique) time i at which 
p*(t) = p and to invert the implicit algebraic equation 
for /bm(P> 0) obtained by setting t = t in Eq. (24). 

Another approach is to evolve the system back- 
wards in time during the time period i defined above, 
from the initial volume fraction _f bm * defining the 
endocortical wall. Indeed, this backwards evolution 
directly provides the bone volume fraction /bm(P>0)- 
This approach is easier to implement numerically as it 



does not require the knowledge of the roots / ; of the 
phenomenological function s v (f). 
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